#!/usr/bin/env python

from scipy import comb, factorial

def binomialCoefficient(n, k):
    if k < 0 or k > n:
        return 0
    if k > n - k: # take advantage of symmetry
        k = n - k
    c = 1
    for i in range(k):
        c = c * (n - (k - (i+1)))
        c = c // (i+1)
    return c


def basis_fun_eval(j,n,a,b,x):
#    if(j > n or j < 0):
#        return 0. # important for derivative evaluation
#    else:
        return binomialCoefficient(n,j)*(x-a)**j*(b-x)**(n-j)/(b-a)**n
#        return comb(n,j,exact=0)*(x-a)**j*(b-x)**(n-j)/(b-a)**n



def basis_fun_derivative(k,j,n,a,b,x):
	if (k > n): 
		print 'bernstein_basis_fun_derivative::Error: k > n!'
	elif (k == 1):
		return n/(b-a)*( basis_fun_eval(j-1,n-1,a,b,x) - basis_fun_eval(j,n-1,a,b,x) ) 
	else:
		return n/(b-a)*(basis_fun_derivative(k-1,j-1,n-1,a,b,x)-basis_fun_derivative(k-1,j,n-1,a,b,x))

def basis_fun_der(p,i,n,a,b,x):
	ks = max(0,i+p-n)
        ke = min(i,p)
        bmar = 1/(b-a)**p
	suma = 0.
        for k in range(ks,ke+1):
#		suma += (-1)**(k+p)*comb(p,k,exact=0)*basis_fun_eval(i-k,n-p,a,b,x)
		suma += (-1)**(k+p)*binomialCoefficient(p,k)*basis_fun_eval(i-k,n-p,a,b,x)
	return float(factorial(n,exact=1))/float(factorial(n-p,exact=1))*bmar*suma

def interpolant(beta,n,a,b,x):  
	f = 0.
	for j in range(n+1):
		f += beta[j]*basis_fun_eval(j,n,a,b,x)
	return f

### Routines for two-dimensional problems

def interpolant2D(n,m,x1,x2,y1,y2,beta,x,y):
    f = 0.
    for i in range(n+1):
        for j in range(m+1):
            f += beta[i][j]*basis_fun_eval(j,m,y1,y2,y)*basis_fun_eval(i,n,x1,x2,x)
    return f

def laplacian(i,j,n,m,x1,x2,y1,y2,x,y):
# Laplacian operator:
    return basis_fun_der(2,i,n,x1,x2,x)*basis_fun_eval(j,m,y1,y2,y) + \
           basis_fun_eval(i,n,x1,x2,x)*basis_fun_der(2,j,m,y1,y2,y)

def lapl(i,j,n,m,x1,x2,y1,y2,x,y):
# Laplacian operator:

    term1 = 0.
    for i in range(n+1):
        for j in range(m+1):
            term1 += basis_fun_derivative(2,i,n,x1,x2,x)*basis_fun_eval(j,m,y1,y2,y)
    term2 = 0.
    for i in range(n+1):
        for j in range(m+1):
            term2 +=basis_fun_eval(i,n,x1,x2,x)*basis_fun_derivative(2,j,m,y1,y2,y)

    return term1+term2
